Investigating Cambodia’s Nightlights Variability

library(Rnightlights)
library(lubridate)
library(reshape2)
library(ggplot2)

pkgOptions(downloadMethod="aria")
pkgOptions(extractMethod = "rast", numCores=4)
pkgOptions(deleteTiles = TRUE)
ctry <- "KHM"

#Getting the nightlights data
#2012
df_Country12 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201301"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces12 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201301"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts12 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201301"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions12 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201301"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages12 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201204"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2013
#the feb tile seems to be corrupted
df_Country13 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201312"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces13 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201312"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts13 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201312"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions13 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201312"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages13 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201312"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2014 
df_Country14 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201412"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces14 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201412"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts14 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201412"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions14 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201412"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages14 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201412"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2015
df_Country15 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201501", "201512"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces15 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201501", "201512"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts15 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201501", "201512"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions15 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201501", "201512"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages15 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201501", "201512"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2016
df_Country16 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201601", "201612"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces16 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201601", "201612"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts16 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201601", "201612"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions16 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201601", "201612"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages16 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201601", "201612"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2017
df_Country17 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201701", "201712"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces17 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201701", "201712"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts17 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201701", "201712"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions17 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201701", "201712"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages17 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201701", "201712"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#2018
#June dataset is missing
df_Country18 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201801", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces18 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201801", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts18 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201801", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions18 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201801", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages18 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201801", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

#July 2018 to April 2019
df_Country19 <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201807", "201904"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Provinces19 <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201807", "201904"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Districts19 <- getCtryNlData(ctryCode = ctry, admLevel = "admin2", nlType = "VIIRS.M", nlPeriods = nlRange("201807", "201904"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Divisions19 <- getCtryNlData(ctryCode = ctry, admLevel = "admin3", nlType = "VIIRS.M", nlPeriods = nlRange("201807", "201904"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))
df_Villages19 <- getCtryNlData(ctryCode = ctry, admLevel = "admin4", nlType = "VIIRS.M", nlPeriods = nlRange("201807", "201904"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

## Melting the data per year
#Its time to get a melt function to avoid copy pasting code
meltData <- function(x){
  x <- melt(x,
            id.vars = grep("NL_", names(x), 
                           invert=TRUE), 
            variable.name = "nlPeriod", 
            value.name = "radiancesum")
  #extract date from the NL col names
  x$nlPeriod <- substr(x$nlPeriod, 12, 17)
  
  #format period as date
  x$nlPeriod <- ymd(paste0(substr(x$nlPeriod, 1,4), 
                           "-",substr(x$nlPeriod, 5,6), "-01"))
  return(x)
}

#melt the data at country level
KHM <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201904"), ignoreMissing=TRUE, nlStats = list("sum",na.rm=TRUE))

KHM <- meltData(KHM)
#time series at country level
ggplot(data = KHM, 
       aes(x=nlPeriod, y=radiancesum, 
           color=KHM[[3]])) +
  scale_x_date(date_breaks = "6 months", date_labels = "%Y-%m")+
  geom_line(color = "#00AFBB", size = 1) + geom_point() + labs(color = names(KHM)[3]) + 
  xlab("Month") + 
  ylab("Sum of Radiances") +
  ggtitle(paste0(unique(names(KHM)[3]), " sum of radiances for ", ctry))

#melt the data at province level
khm12 <- meltData(df_Provinces12)
khm13 <- meltData(df_Provinces13)
khm14 <- meltData(df_Provinces14)
khm15 <- meltData(df_Provinces15)
khm16 <- meltData(df_Provinces16)
khm17 <- meltData(df_Provinces17)
khm18 <- meltData(df_Provinces18)
khm19 <- meltData(df_Provinces19)

#time series at province level
KHMProv <- getCtryNlData(ctryCode = ctry, admLevel = "admin1", nlType = "VIIRS.M", nlPeriods = nlRange("201204", "201904"), ignoreMissing=TRUE, nlStats = list("sum",na.rm=TRUE))
KHMProv <- meltData(KHMProv)

ggplot(data = KHMProv, 
       aes(x=nlPeriod, y=radiancesum, 
           color=KHMProv[[2]])) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b-%Y")+
  geom_line() + geom_point() + labs(color = names(KHMProv)[2]) + 
  xlab("Month") + 
  ylab("Sum of Radiances") +
  ggtitle(paste0(unique(names(KHMProv)[2]), " sum of radiances for ", ctry))

#ts' plot is hard to discern, let's facet_wrap & see yearly trends per province
level112 <- ggplot(data = khm12) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 ) 
level112

level113 <- ggplot(data = khm13) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 ) 
level113

level114 <- ggplot(data = khm14) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level114

level115 <- ggplot(data = khm15) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level115

level116 <- ggplot(data = khm16) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level116

level117 <- ggplot(data = khm17) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level117

level118 <- ggplot(data = khm18) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level118

level119 <- ggplot(data = khm19) +
  geom_line(mapping = aes(x = nlPeriod, y = radiancesum)) +
  facet_wrap(~ province_.khêt., nrow = 5 )
level119

#2013 had low nl but we don't have the feb data
#See article:hot dark season of 2013
#http://www.sebastianstrangio.com/2014/05/27/its-electricity-dwindling-cambodia-is-getting-very-dark-and-very-very-hot/
 
#data wrangling time
#at village level
KHMVil <- getCtryNlData(ctryCode = ctry, admLevel = "admin0", nlType = "VIIRS.M", nlPeriods = nlRange("201303", "201805"), ignoreMissing=FALSE, nlStats = list("sum",na.rm=TRUE))

KHMVil <- meltData(KHMVil)

#2015, notice the way values dip around July?
#investigate at village level, which villages have 0 values
#Does the lower nl provinces contribute to the zeros??Not really
#Does normalization cause any difference??

library(dplyr)
# names(KHMVil)
#  bb <- KHMVil %>%
#    group_by(province_.khêt.) %>%
#    filter(radiancesum <= 1) %>%
#    summarise(n = n()) %>%
#    arrange(desc(n))